import java.util.Random;
import java.util.Arrays;
import java.lang.Math;
import java.io.*;

/*
Comparing error (absolute) of each sketch as memory percent changes
Varying size of (both) data sets
Varying size of one data set (keeping other fixed) to measure skew
Varying the actual difference between streams (scatter plot)
*/

public class TwoSampleKStest
{
    public static void main(String[] args) throws Exception
    {
	if (args.length < 2)
	{
	    System.out.println("Insufficient arguments");
	    System.exit(0);
	}

	TwoSampleKStest ks = new TwoSampleKStest();

	if (args[0].equalsIgnoreCase("M"))  // vary memory
	{	
	    if (args.length > 3)
		ks.varyMemory(args[1], args[2], args[3]);
	    else
		ks.varyMemory(args[1], args[2], "results/temp-" + args[0].substring(args[0].lastIndexOf("/") + 1) + args[1].substring(args[1].lastIndexOf("/") + 1) + "results");
	}
	else if (args[0].equalsIgnoreCase("B"))  // vary both sizes
	{
	    if (args.length > 3)
                ks.varyBothdatasize(args[1], args[2], args[3]);
            else
                ks.varyBothdatasize(args[1], args[2], "results/temp-" + args[0].substring(args[0].lastIndexOf("/") + 1) + args[1].substring(args[1].lastIndexOf("/") + 1) + "results");
	}
	else if (args[0].equalsIgnoreCase("O"))  // vary both sizes
	{
	    if (args.length > 3)
		ks.varyOnedatasize(args[1], args[2], args[3]);
	    else
		ks.varyOnedatasize(args[1], args[2], "results/temp-" + args[0].substring(args[0].lastIndexOf("/") + 1) + args[1].substring(args[1].lastIndexOf("/") + 1) + "results");
	}
	else if (args[0].equalsIgnoreCase("D"))  // collect real/est pairs for many different distances
	{   
	    if (args.length > 3)
		ks.varyDistance(args[1], args[2], args[3]);
	    else
		ks.varyDistance(args[1], args[2], "results/temp-" + args[0].substring(args[0].lastIndexOf("/") + 1) + args[1].substring(args[1].lastIndexOf("/") + 1) + "results");
	}
	else if (args[0].equalsIgnoreCase("R"))
	{   
	    if (args.length > 3)
		ks.varyMemoryReal(args[1], args[2], args[3]);
	    else
		ks.varyMemoryReal(args[1], args[2], "results/temp-" + args[0].substring(args[0].lastIndexOf("/") + 1) + args[1].substring(args[1].lastIndexOf("/") + 1) + "results");
	}

    }

    // Comparing error (absolute) of each sketch as memory percent changes 
    public void varyMemory(String d1, String d2, String results) throws Exception // to simplify file reading and number parsing
    {
	int n = 10000, m = 10000;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	double[] frac = {0.001, 0.0031622, 0.01, 0.031622, 0.1}; //{0.05, 0.1, 0.15, 0.20, 0.25};
	int fracs = frac.length;
	double[][] relErr = new double[4][fracs];
	double[][] absErr = new double[4][fracs];
	
	for (int i = 0; i < fracs; ++i)
	    for (int j = 0; j < 4; ++j)
		relErr[j][i] = absErr[j][i] = 0.0; 
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int k = 0; k < datarepeats; ++k)
	{
	    // read data from file
	    BufferedReader r1 = new BufferedReader(new FileReader(d1 + k + ".dat"));
	    n = Integer.parseInt(r1.readLine());
	    double[] data1 = new double[n];
	    for (int i = 0; i < n; i++)
		data1[i] = Double.parseDouble(r1.readLine());
	    r1.close();
	    
	    BufferedReader r2 = new BufferedReader(new FileReader(d2 + k + ".dat"));
	    m = Integer.parseInt(r2.readLine());
	    double[] data2 = new double[m];
	    for (int i = 0; i < m; i++)
		data2[i] = Double.parseDouble(r2.readLine());
	    r2.close();
	    
	    // compute the real KS statistic
	    double D_real = exactKS(data1, data2); 
	    System.out.println("Real D: " + D_real);
	    w.write("Real D: " + D_real + "\n");
	    	    
	    for (int l = 0; l < fracs; ++l)
	    {
		double f = frac[l];
		System.out.println("\nf = " + f);
		w.write("\nf = " + f + "\n");

		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary gk1 = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    gk1.insert(data1[i]);
		
		GreenwaldKhannaQuantileSummary gk2 = new GreenwaldKhannaQuantileSummary((int)(f * m));
		for (int i = 0; i < m; ++i)
		    gk2.insert(data2[i]);
		
		// compute the estimated KS statistic from GK sketch
		double D_GK = estimateKS(gk1, n, gk2, m);

		System.out.print("GK: " + D_GK);
		w.write("GK: " + D_GK);
		System.out.print(", rel error: " + relError(D_GK, D_real));
		w.write(", rel error: " + relError(D_GK, D_real));
		System.out.println(", abs error: " + Math.abs(D_GK - D_real));
		w.write(", abs error: " + Math.abs(D_GK - D_real) + "\n");
		relErr[0][l] += relError(D_GK, D_real)/datarepeats;
		absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
	    
	    
		/** test the q-Digest data structure **/
		QDigest qd1 = new QDigest((int)(f * n)/3);   // compression factor = memory/3
		for (int i = 0; i < n; ++i)
		    qd1.insert(data1[i]);
		
		QDigest qd2 = new QDigest((int)(f * m)/3);
		for (int i = 0; i < m; ++i)
		    qd2.insert(data2[i]);

		// compute the estimated KS statistic from QDigest sketch
		double D_QD = estimateKS(qd1, n, qd2, m);

                System.out.print("QD: " + D_QD);
                w.write("QD: " + D_QD);
		System.out.print(", rel error: " + relError(D_QD, D_real));
                w.write(", rel error: " + relError(D_QD, D_real));
		System.out.println(", abs error: " + Math.abs(D_QD - D_real));
		w.write(", abs error: " + Math.abs(D_QD - D_real) + "\n");
		relErr[1][l] += relError(D_QD, D_real)/datarepeats;
                absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats;
	    

		/** test the reservoir sampling data structure **/
		ReservoirSampler rs1 = new ReservoirSampler((int)(f * n));
		for (int i = 0; i < n; ++i)
		    rs1.insert(data1[i]);
		
		ReservoirSampler rs2 = new ReservoirSampler((int)(f * m));
		for (int i = 0; i < m; ++i)
		    rs2.insert(data2[i]);

		// compute the estimated KS statistic from reservoir sampler
		double D_RS = estimateKS(rs1, n, rs2, m);

                System.out.print("RS: " + D_RS);
		w.write("RS: " + D_RS);
                System.out.print(", rel error: " + relError(D_RS, D_real));
                w.write(", rel error: " + relError(D_RS, D_real));
		System.out.println(", abs error: " + Math.abs(D_RS - D_real));
		w.write(", abs error: " + Math.abs(D_RS - D_real) + "\n");
		relErr[2][l] += relError(D_RS, D_real)/datarepeats;
                absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats;

		/** Compute value from exact KS on sampled data **/
		double[] sample1 = rs1.getSamples();
		double[] sample2 = rs2.getSamples();
		double D_SA = exactKS(sample1, sample2);

		System.out.print("SA: " + D_SA);
                w.write("SA: " + D_SA);
                System.out.print(", rel error: " + relError(D_SA, D_real));
                w.write(", rel error: " + relError(D_SA, D_real));
                System.out.println(", abs error: " + Math.abs(D_SA - D_real));
                w.write(", abs error: " + Math.abs(D_SA - D_real) + "\n");
                relErr[3][l] += relError(D_SA, D_real)/datarepeats;
                absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    }

	    w.flush();
	}

	System.out.println("\n\n");
	for (int i = 0; i < fracs; ++i)
	{
	    System.out.println("f = " + frac[i]);
	    System.out.println("GK rel: " + relErr[0][i]);
	    System.out.println("GK abs: " + absErr[0][i]);
	    System.out.println("QD rel: " + relErr[1][i]);
            System.out.println("QD abs: " + absErr[1][i]);
	    System.out.println("RS rel: " + relErr[2][i]);
            System.out.println("RS abs: " + absErr[2][i]);
	    System.out.println("SA rel: " + relErr[3][i]);
            System.out.println("SA abs: " + absErr[3][i]);
	    System.out.println("");

	    w.write("f = " + frac[i] + "\n");
            w.write("GK rel: " + relErr[0][i] + "\n");
            w.write("GK abs: " + absErr[0][i] + "\n");
            w.write("QD rel: " + relErr[1][i] + "\n");
            w.write("QD abs: " + absErr[1][i] + "\n");
            w.write("RS rel: " + relErr[2][i] + "\n");
            w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
            w.write("SA abs: " + absErr[3][i] + "\n");
            w.write("\n");

	}
	w.close();

	w = new BufferedWriter(new FileWriter(results + ".dat"));
	for (int i = 0; i < fracs; ++i)
	{
	    w.write(frac[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
	w.close();
    }

    // Varying size of (both) data sets
    public void varyBothdatasize(String d1, String d2, String results) throws Exception // to simplify file reading and number parsing
    {
	int[] ens = {1000, 10000, 100000};
	int num = ens.length;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	double f = 0.01;
	double[][] relErr = new double[4][num];
	double[][] absErr = new double[4][num];
	
	for (int i = 0; i < num; ++i)
	    for (int j = 0; j < 4; ++j)
		relErr[j][i] = absErr[j][i] = 0.0; 
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int l = 0; l < ens.length; ++l)
	{
	    int n = ens[l], m = n;

	    for (int k = 0; k < datarepeats; ++k)
	    {
		// read data from file
		BufferedReader r1 = new BufferedReader(new FileReader(d1 + n + "-" + k + ".dat"));
		if (Integer.parseInt(r1.readLine()) != n) { System.out.println("data mismatch!"); System.exit(0); }
		double[] data1 = new double[n];
		for (int i = 0; i < n; i++)
		    data1[i] = Double.parseDouble(r1.readLine());
		r1.close();
	    
		BufferedReader r2 = new BufferedReader(new FileReader(d2 + m + "-" + k + ".dat"));
		if (Integer.parseInt(r2.readLine()) != m) { System.out.println("data mismatch!"); System.exit(0); }
		double[] data2 = new double[m];
		for (int i = 0; i < m; i++)
		    data2[i] = Double.parseDouble(r2.readLine());
		r2.close();
	    
		// compute the real KS statistic
		double D_real = exactKS(data1, data2); 
		System.out.println("Real D: " + D_real);
		w.write("Real D: " + D_real + "\n");

		System.out.println("\nn = " + n);
		w.write("\nn = " + n + "\n");

		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary gk1 = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    gk1.insert(data1[i]);
		
		GreenwaldKhannaQuantileSummary gk2 = new GreenwaldKhannaQuantileSummary((int)(f * m));
		for (int i = 0; i < m; ++i)
		    gk2.insert(data2[i]);
		
		// compute the estimated KS statistic from GK sketch
		double D_GK = estimateKS(gk1, n, gk2, m);

		System.out.print("GK: " + D_GK);
		w.write("GK: " + D_GK);
		System.out.print(", rel error: " + relError(D_GK, D_real));
		w.write(", rel error: " + relError(D_GK, D_real));
		System.out.println(", abs error: " + Math.abs(D_GK - D_real));
		w.write(", abs error: " + Math.abs(D_GK - D_real) + "\n");
		relErr[0][l] += relError(D_GK, D_real)/datarepeats;
		absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
	    
	    
		/** test the q-Digest data structure **/
		QDigest qd1 = new QDigest((int)(f * n)/3);   // compression factor = memory/3
		for (int i = 0; i < n; ++i)
		    qd1.insert(data1[i]);
		
		QDigest qd2 = new QDigest((int)(f * m)/3);
		for (int i = 0; i < m; ++i)
		    qd2.insert(data2[i]);
		
		// compute the estimated KS statistic from QDigest sketch
		double D_QD = estimateKS(qd1, n, qd2, m);

                System.out.print("QD: " + D_QD);
                w.write("QD: " + D_QD);
		System.out.print(", rel error: " + relError(D_QD, D_real));
                w.write(", rel error: " + relError(D_QD, D_real));
		System.out.println(", abs error: " + Math.abs(D_QD - D_real));
		w.write(", abs error: " + Math.abs(D_QD - D_real) + "\n");
		relErr[1][l] += relError(D_QD, D_real)/datarepeats;
                absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats;
	    

		/** test the reservoir sampling data structure **/
		ReservoirSampler rs1 = new ReservoirSampler((int)(f * n));
		for (int i = 0; i < n; ++i)
		    rs1.insert(data1[i]);
		
		ReservoirSampler rs2 = new ReservoirSampler((int)(f * m));
		for (int i = 0; i < m; ++i)
		    rs2.insert(data2[i]);

		// compute the estimated KS statistic from reservoir sampler
		double D_RS = estimateKS(rs1, n, rs2, m);

                System.out.print("RS: " + D_RS);
		w.write("RS: " + D_RS);
                System.out.print(", rel error: " + relError(D_RS, D_real));
                w.write(", rel error: " + relError(D_RS, D_real));
		System.out.println(", abs error: " + Math.abs(D_RS - D_real));
		w.write(", abs error: " + Math.abs(D_RS - D_real) + "\n");
		relErr[2][l] += relError(D_RS, D_real)/datarepeats;
                absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats;

		/** Compute value from exact KS on sampled data **/
		double[] sample1 = rs1.getSamples();
		double[] sample2 = rs2.getSamples();
		double D_SA = exactKS(sample1, sample2);

		System.out.print("SA: " + D_SA);
                w.write("SA: " + D_SA);
                System.out.print(", rel error: " + relError(D_SA, D_real));
                w.write(", rel error: " + relError(D_SA, D_real));
                System.out.println(", abs error: " + Math.abs(D_SA - D_real));
                w.write(", abs error: " + Math.abs(D_SA - D_real) + "\n");
                relErr[3][l] += relError(D_SA, D_real)/datarepeats;
                absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    }

	    w.flush();
	}

	System.out.println("\n\n");
	for (int i = 0; i < num; ++i)
	{
	    System.out.println("n = " + ens[i]);
	    System.out.println("GK rel: " + relErr[0][i]);
	    System.out.println("GK abs: " + absErr[0][i]);
	    System.out.println("QD rel: " + relErr[1][i]);
            System.out.println("QD abs: " + absErr[1][i]);
	    System.out.println("RS rel: " + relErr[2][i]);
            System.out.println("RS abs: " + absErr[2][i]);
	    System.out.println("SA rel: " + relErr[3][i]);
            System.out.println("SA abs: " + absErr[3][i]);
	    System.out.println("");

	    w.write("n = " + ens[i] + "\n");
            w.write("GK rel: " + relErr[0][i] + "\n");
            w.write("GK abs: " + absErr[0][i] + "\n");
            w.write("QD rel: " + relErr[1][i] + "\n");
            w.write("QD abs: " + absErr[1][i] + "\n");
            w.write("RS rel: " + relErr[2][i] + "\n");
            w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
            w.write("SA abs: " + absErr[3][i] + "\n");
            w.write("\n");

	}

	w.close();

	w = new BufferedWriter(new FileWriter(results + ".dat"));
        for (int i = 0; i < num; ++i)
	{
	    w.write(ens[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
        w.close();

    }

    public void varyOnedatasize(String d1, String d2, String results) throws Exception // to simplify file reading and number parsing
    {
	int n = 10000;
	int[] ens = {1000, 10000, 100000};
	int num = ens.length;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	double f = 0.01;
	double[][] relErr = new double[4][num];
	double[][] absErr = new double[4][num];
	
	for (int i = 0; i < num; ++i)
	    for (int j = 0; j < 4; ++j)
		relErr[j][i] = absErr[j][i] = 0.0; 
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int l = 0; l < ens.length; ++l)
	{
	    int m = ens[l];

	    for (int k = 0; k < datarepeats; ++k)
	    {
		// read data from file
		BufferedReader r1 = new BufferedReader(new FileReader(d1 + n + "-" + k + ".dat"));
		if (Integer.parseInt(r1.readLine()) != n) { System.out.println("data mismatch!"); System.exit(0); }
		double[] data1 = new double[n];
		for (int i = 0; i < n; i++)
		    data1[i] = Double.parseDouble(r1.readLine());
		r1.close();
	    
		BufferedReader r2 = new BufferedReader(new FileReader(d2 + m + "-" + k + ".dat"));
		if (Integer.parseInt(r2.readLine()) != m) { System.out.println("data mismatch!"); System.exit(0); }
		double[] data2 = new double[m];
		for (int i = 0; i < m; i++)
		    data2[i] = Double.parseDouble(r2.readLine());
		r2.close();
	    
		// compute the real KS statistic
		double D_real = exactKS(data1, data2); 
		System.out.println("Real D: " + D_real);
		w.write("Real D: " + D_real + "\n");

		System.out.println("\nn = " + n + " m = " + m);
		w.write("\nn = " + n + " m = " + m + "\n");

		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary gk1 = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    gk1.insert(data1[i]);
		
		GreenwaldKhannaQuantileSummary gk2 = new GreenwaldKhannaQuantileSummary((int)(f * m));
		for (int i = 0; i < m; ++i)
		    gk2.insert(data2[i]);
		
		// compute the estimated KS statistic from GK sketch
		double D_GK = estimateKS(gk1, n, gk2, m);

		System.out.print("GK: " + D_GK);
		w.write("GK: " + D_GK);
		System.out.print(", rel error: " + relError(D_GK, D_real));
		w.write(", rel error: " + relError(D_GK, D_real));
		System.out.println(", abs error: " + Math.abs(D_GK - D_real));
		w.write(", abs error: " + Math.abs(D_GK - D_real) + "\n");
		relErr[0][l] += relError(D_GK, D_real)/datarepeats;
		absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
	    
	    
		/** test the q-Digest data structure **/
		QDigest qd1 = new QDigest((int)(f * n)/3);   // compression factor = memory/3
		for (int i = 0; i < n; ++i)
		    qd1.insert(data1[i]);
		
		QDigest qd2 = new QDigest((int)(f * m)/3);
		for (int i = 0; i < m; ++i)
		    qd2.insert(data2[i]);
		
		// compute the estimated KS statistic from QDigest sketch
		double D_QD = estimateKS(qd1, n, qd2, m);

                System.out.print("QD: " + D_QD);
                w.write("QD: " + D_QD);
		System.out.print(", rel error: " + relError(D_QD, D_real));
                w.write(", rel error: " + relError(D_QD, D_real));
		System.out.println(", abs error: " + Math.abs(D_QD - D_real));
		w.write(", abs error: " + Math.abs(D_QD - D_real) + "\n");
		relErr[1][l] += relError(D_QD, D_real)/datarepeats;
                absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats;
	    

		/** test the reservoir sampling data structure **/
		ReservoirSampler rs1 = new ReservoirSampler((int)(f * n));
		for (int i = 0; i < n; ++i)
		    rs1.insert(data1[i]);
		
		ReservoirSampler rs2 = new ReservoirSampler((int)(f * m));
		for (int i = 0; i < m; ++i)
		    rs2.insert(data2[i]);

		// compute the estimated KS statistic from reservoir sampler
		double D_RS = estimateKS(rs1, n, rs2, m);

                System.out.print("RS: " + D_RS);
		w.write("RS: " + D_RS);
                System.out.print(", rel error: " + relError(D_RS, D_real));
                w.write(", rel error: " + relError(D_RS, D_real));
		System.out.println(", abs error: " + Math.abs(D_RS - D_real));
		w.write(", abs error: " + Math.abs(D_RS - D_real) + "\n");
		relErr[2][l] += relError(D_RS, D_real)/datarepeats;
                absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats;

		/** Compute value from exact KS on sampled data **/
		double[] sample1 = rs1.getSamples();
		double[] sample2 = rs2.getSamples();
		double D_SA = exactKS(sample1, sample2);

		System.out.print("SA: " + D_SA);
                w.write("SA: " + D_SA);
                System.out.print(", rel error: " + relError(D_SA, D_real));
                w.write(", rel error: " + relError(D_SA, D_real));
                System.out.println(", abs error: " + Math.abs(D_SA - D_real));
                w.write(", abs error: " + Math.abs(D_SA - D_real) + "\n");
                relErr[3][l] += relError(D_SA, D_real)/datarepeats;
                absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    }

	    w.flush();
	}

	System.out.println("\n\n");
	for (int i = 0; i < num; ++i)
	{
	    System.out.println("n = " + ens[i]);
	    System.out.println("GK rel: " + relErr[0][i]);
	    System.out.println("GK abs: " + absErr[0][i]);
	    System.out.println("QD rel: " + relErr[1][i]);
            System.out.println("QD abs: " + absErr[1][i]);
	    System.out.println("RS rel: " + relErr[2][i]);
            System.out.println("RS abs: " + absErr[2][i]);
	    System.out.println("SA rel: " + relErr[3][i]);
            System.out.println("SA abs: " + absErr[3][i]);
	    System.out.println("");

	    w.write("n = " + ens[i] + "\n");
            w.write("GK rel: " + relErr[0][i] + "\n");
            w.write("GK abs: " + absErr[0][i] + "\n");
            w.write("QD rel: " + relErr[1][i] + "\n");
            w.write("QD abs: " + absErr[1][i] + "\n");
            w.write("RS rel: " + relErr[2][i] + "\n");
            w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
            w.write("SA abs: " + absErr[3][i] + "\n");
            w.write("\n");

	}

	w.close();

	w = new BufferedWriter(new FileWriter(results + ".dat"));
        for (int i = 0; i < num; ++i)
	{
	    w.write(ens[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
        w.close();

    }

    public void varyDistance(String d1, String d2, String results) throws Exception // to simplify file reading and number parsing
    {
	int n = 10000, m = 10000;
	double f = 0.01;
	Random r = new Random();
	String[] mean = {"0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1"}; // used strings to preserve formatting 
	int means = mean.length;
	double[] real = new double[means];
	double[][] est = new double[4][means];
	int c = 0;

	for (int i = 0; i < means; ++i)
	{
	    real[i] = 0.0;
	    for (int j = 0; j < 4; ++j)
		est[j][i] = 0.0;
	}
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	// read data from file
	BufferedReader r1 = new BufferedReader(new FileReader(d1));
	n = Integer.parseInt(r1.readLine());
	double[] data1 = new double[n];
	for (int i = 0; i < n; i++)
	    data1[i] = Double.parseDouble(r1.readLine());
	r1.close();
	 
	for (int l = 0; l < means; ++l)
	{
	    System.out.println("mean = " + mean[l]);
	    w.write("\nmean = " + mean[l] + "\n"); 

	    BufferedReader r2 = new BufferedReader(new FileReader(d2 + mean[l] +  "-1.dat"));
	    m = Integer.parseInt(r2.readLine());
	    double[] data2 = new double[m];
	    for (int i = 0; i < m; i++)
		data2[i] = Double.parseDouble(r2.readLine());
	    r2.close();
	    
	    // compute the real KS statistic
	    double D_real = exactKS(data1, data2); 
	    System.out.println("Real D: " + D_real);
	    w.write("Real D: " + D_real + "\n");
	    real[c] = D_real;
	    
	    System.out.println("\nf = " + f);
	    w.write("\nf = " + f + "\n");
	    
	    /** test the GK-data structure **/
	    GreenwaldKhannaQuantileSummary gk1 = new GreenwaldKhannaQuantileSummary((int)(f * n));
	    for (int i = 0; i < n; ++i)
		gk1.insert(data1[i]);
	    
	    GreenwaldKhannaQuantileSummary gk2 = new GreenwaldKhannaQuantileSummary((int)(f * m));
	    for (int i = 0; i < m; ++i)
		gk2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from GK sketch
	    double D_GK = estimateKS(gk1, n, gk2, m);
	    
	    System.out.print("GK: " + D_GK);
	    w.write("GK: " + D_GK);
	    System.out.print(", rel error: " + relError(D_GK, D_real));
	    w.write(", rel error: " + relError(D_GK, D_real));
	    System.out.println(", abs error: " + Math.abs(D_GK - D_real));
	    w.write(", abs error: " + Math.abs(D_GK - D_real) + "\n");
	    est[0][c] = D_GK;
	    
	    
	    /** test the q-Digest data structure **/
	    QDigest qd1 = new QDigest((int)(f * n)/3);   // compression factor = memory/3
	    for (int i = 0; i < n; ++i)
		qd1.insert(data1[i]);
	    
	    QDigest qd2 = new QDigest((int)(f * m)/3);
	    for (int i = 0; i < m; ++i)
		qd2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from QDigest sketch
	    double D_QD = estimateKS(qd1, n, qd2, m);
	    
	    System.out.print("QD: " + D_QD);
	    w.write("QD: " + D_QD);
	    System.out.print(", rel error: " + relError(D_QD, D_real));
	    w.write(", rel error: " + relError(D_QD, D_real));
	    System.out.println(", abs error: " + Math.abs(D_QD - D_real));
	    w.write(", abs error: " + Math.abs(D_QD - D_real) + "\n");
	    est[1][c] = D_QD;
	    
	    /** test the reservoir sampling data structure **/
	    ReservoirSampler rs1 = new ReservoirSampler((int)(f * n));
	    for (int i = 0; i < n; ++i)
		rs1.insert(data1[i]);
	    
	    ReservoirSampler rs2 = new ReservoirSampler((int)(f * m));
	    for (int i = 0; i < m; ++i)
		rs2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from reservoir sampler
	    double D_RS = estimateKS(rs1, n, rs2, m);
	    
	    System.out.print("RS: " + D_RS);
	    w.write("RS: " + D_RS);
	    System.out.print(", rel error: " + relError(D_RS, D_real));
	    w.write(", rel error: " + relError(D_RS, D_real));
	    System.out.println(", abs error: " + Math.abs(D_RS - D_real));
	    w.write(", abs error: " + Math.abs(D_RS - D_real) + "\n");
	    est[2][c] = D_RS;
	    
	    /** Compute value from exact KS on sampled data **/
	    double[] sample1 = rs1.getSamples();
	    double[] sample2 = rs2.getSamples();
	    double D_SA = exactKS(sample1, sample2);
	    
	    System.out.print("SA: " + D_SA);
	    w.write("SA: " + D_SA);
	    System.out.print(", rel error: " + relError(D_SA, D_real));
	    w.write(", rel error: " + relError(D_SA, D_real));
	    System.out.println(", abs error: " + Math.abs(D_SA - D_real));
	    w.write(", abs error: " + Math.abs(D_SA - D_real) + "\n");
	    est[3][c] = D_SA;
	    
	    c = c + 1;
	    w.flush();
	}

	System.out.println("\n\n");
	for (int i = 0; i < means; ++i)
        {
	    System.out.print(real[i]);
	    for (int j = 0; j < 4; ++j)
		System.out.print("\t" + est[j][i]);
	    System.out.println("");
	}
	
	w.close();
	
	w = new BufferedWriter(new FileWriter(results + ".dat"));
	for (int i = 0; i < means; ++i)
	{
	    w.write(real[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + est[j][i]);
	    w.write("\n");
	}
	w.close();
    }

   // Comparing error (absolute) of each sketch on REAL data as memory percent changes 
    public void varyMemoryReal(String d1, String d2, String results) throws Exception // to simplify file reading and number parsing
    {
	int n = 10000, m = 10000;
	int repeats = 1, datarepeats = 1;
	Random r = new Random();
	double[] frac = {0.001, 0.0031622, 0.01, 0.031622, 0.1}; //{0.05, 0.1, 0.15, 0.20, 0.25};
	int fracs = frac.length;
	double[][] relErr = new double[4][fracs];
	double[][] absErr = new double[4][fracs];
	
	for (int i = 0; i < fracs; ++i)
	    for (int j = 0; j < 4; ++j)
		relErr[j][i] = absErr[j][i] = 0.0; 
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	// read data from file
	BufferedReader r1 = new BufferedReader(new FileReader(d1));
	n = Integer.parseInt(r1.readLine());
	double[] data1 = new double[n];
	for (int i = 0; i < n; i++)
	    data1[i] = Double.parseDouble(r1.readLine());
	r1.close();
	
	BufferedReader r2 = new BufferedReader(new FileReader(d2));
	m = Integer.parseInt(r2.readLine());
	double[] data2 = new double[m];
	for (int i = 0; i < m; i++)
	    data2[i] = Double.parseDouble(r2.readLine());
	r2.close();
	
	// compute the real KS statistic
	double D_real = exactKS(data1, data2); 
	System.out.println("Real D: " + D_real);
	w.write("Real D: " + D_real + "\n");
	
	for (int l = 0; l < fracs; ++l)
	{
	    double f = frac[l];
	    System.out.println("\nf = " + f);
	    w.write("\nf = " + f + "\n");
	    
	    /** test the GK-data structure **/
	    GreenwaldKhannaQuantileSummary gk1 = new GreenwaldKhannaQuantileSummary((int)(f * n));
	    for (int i = 0; i < n; ++i)
		gk1.insert(data1[i]);
	    
	    GreenwaldKhannaQuantileSummary gk2 = new GreenwaldKhannaQuantileSummary((int)(f * m));
	    for (int i = 0; i < m; ++i)
		gk2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from GK sketch
	    double D_GK = estimateKS(gk1, n, gk2, m);
	    
	    System.out.print("GK: " + D_GK);
	    w.write("GK: " + D_GK);
	    System.out.print(", rel error: " + relError(D_GK, D_real));
	    w.write(", rel error: " + relError(D_GK, D_real));
	    System.out.println(", abs error: " + Math.abs(D_GK - D_real));
	    w.write(", abs error: " + Math.abs(D_GK - D_real) + "\n");
	    relErr[0][l] += relError(D_GK, D_real)/datarepeats;
	    absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
	    
	    
	    /** test the q-Digest data structure **/
	    QDigest qd1 = new QDigest((int)(f * n)/3);   // compression factor = memory/3
	    for (int i = 0; i < n; ++i)
		qd1.insert(data1[i]);
	    
	    QDigest qd2 = new QDigest((int)(f * m)/3);
	    for (int i = 0; i < m; ++i)
		qd2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from QDigest sketch
	    double D_QD = estimateKS(qd1, n, qd2, m);
	    
	    System.out.print("QD: " + D_QD);
	    w.write("QD: " + D_QD);
	    System.out.print(", rel error: " + relError(D_QD, D_real));
	    w.write(", rel error: " + relError(D_QD, D_real));
	    System.out.println(", abs error: " + Math.abs(D_QD - D_real));
	    w.write(", abs error: " + Math.abs(D_QD - D_real) + "\n");
	    relErr[1][l] += relError(D_QD, D_real)/datarepeats;
	    absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats;
	    
	    
	    /** test the reservoir sampling data structure **/
	    ReservoirSampler rs1 = new ReservoirSampler((int)(f * n));
	    for (int i = 0; i < n; ++i)
		rs1.insert(data1[i]);
	    
	    ReservoirSampler rs2 = new ReservoirSampler((int)(f * m));
	    for (int i = 0; i < m; ++i)
		rs2.insert(data2[i]);
	    
	    // compute the estimated KS statistic from reservoir sampler
	    double D_RS = estimateKS(rs1, n, rs2, m);
	    
	    System.out.print("RS: " + D_RS);
	    w.write("RS: " + D_RS);
	    System.out.print(", rel error: " + relError(D_RS, D_real));
	    w.write(", rel error: " + relError(D_RS, D_real));
	    System.out.println(", abs error: " + Math.abs(D_RS - D_real));
	    w.write(", abs error: " + Math.abs(D_RS - D_real) + "\n");
	    relErr[2][l] += relError(D_RS, D_real)/datarepeats;
	    absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats;

	    /** Compute value from exact KS on sampled data **/
	    double[] sample1 = rs1.getSamples();
	    double[] sample2 = rs2.getSamples();
	    double D_SA = exactKS(sample1, sample2);
	    
	    System.out.print("SA: " + D_SA);
	    w.write("SA: " + D_SA);
	    System.out.print(", rel error: " + relError(D_SA, D_real));
	    w.write(", rel error: " + relError(D_SA, D_real));
	    System.out.println(", abs error: " + Math.abs(D_SA - D_real));
	    w.write(", abs error: " + Math.abs(D_SA - D_real) + "\n");
	    relErr[3][l] += relError(D_SA, D_real)/datarepeats;
	    absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    
	    w.flush();
	}
	
	System.out.println("\n\n");
	for (int i = 0; i < fracs; ++i)
	{
	    System.out.println("f = " + frac[i]);
	    System.out.println("GK rel: " + relErr[0][i]);
	    System.out.println("GK abs: " + absErr[0][i]);
	    System.out.println("QD rel: " + relErr[1][i]);
            System.out.println("QD abs: " + absErr[1][i]);
	    System.out.println("RS rel: " + relErr[2][i]);
            System.out.println("RS abs: " + absErr[2][i]);
	    System.out.println("SA rel: " + relErr[3][i]);
            System.out.println("SA abs: " + absErr[3][i]);
	    System.out.println("");

	    w.write("f = " + frac[i] + "\n");
            w.write("GK rel: " + relErr[0][i] + "\n");
            w.write("GK abs: " + absErr[0][i] + "\n");
            w.write("QD rel: " + relErr[1][i] + "\n");
            w.write("QD abs: " + absErr[1][i] + "\n");
            w.write("RS rel: " + relErr[2][i] + "\n");
            w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
            w.write("SA abs: " + absErr[3][i] + "\n");
            w.write("\n");

	}
	w.close();

	w = new BufferedWriter(new FileWriter(results + ".dat"));
	for (int i = 0; i < fracs; ++i)
	{
	    w.write(frac[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
	w.close();
    }


    // estimates the KS-distance using the algorithm from the paper
    public static double estimateKS(QuantileSketch a, int n, QuantileSketch b, int m)
    {
	// extract all the elements in both sketches
	double[] X = new double[n + m];
	for (int i = 0; i < n; ++i)
	    X[i] = a.getQuantile(i*1.0/n);
	for (int i = 0; i < m; ++i)
            X[i + n] = b.getQuantile(i*1.0/m);

	double D = 0.0;
	// for each value X, estimate the quantity E_X
	for (int i = 0; i < n + m; i++)
	{
	    int ia = a.reverseQuantile(X[i], n);
	    int ib = b.reverseQuantile(X[i], m);
	    double E_X = Math.abs(((double)ia)/n - ((double)ib)/m);
	    D = Math.max(D, E_X);
	}

	return D;
    }

    // adapted from: http://root.cern.ch/root/html/src/TMath.cxx.html#SX5s7D
    public static double exactKS(double[] d1, double[] d2)
    {
	double D_real = 0.0;
	int n = d1.length, m = d2.length;
	int i1 = 0, i2 = 0;
	double diff = 0.0;
	
	// created sorted version of the data                                                                                             
	double[] data1 = new double[n];
	for (int i = 0; i < n; i++)
	    data1[i] = d1[i];
	Arrays.sort(data1);
	double[] data2 = new double[m];
	for (int i = 0; i < m; i++)
	    data2[i] = d2[i];
	Arrays.sort(data2);

	for (int i = 0; i < n + m; i++)
        {
	    if (data1[i1] < data2[i2])
	    {
		diff -= 1.0/n;
		i1++;
		if (i1 >= n) break;
	    }
	    else if (data1[i1] > data2[i2])
	    {
		diff += 1.0/m;
		i2++;
		if (i2 >= m) break;
	    }
	    else
	    {
		double x = data1[i1];
		while (i1 < n && data1[i1] == x)
		{
		    diff -= 1.0/n;
		    i1++;
		}
		while (i2 < m && data2[i2] == x)
		{
		    diff += 1.0/m;
		    i2++;
		}
		if (i1 >= n || i2 >= m) break;
	    }

	    D_real = Math.max(D_real, Math.abs(diff));
	}

	return D_real;
    }
    

    // computes relative error of est from real
    public static double relError(double est, double real)
    {
	return Math.abs(est - real)/real;
    }
    
    // the cdf of N(mu, sigma) on value x
    public static double cdfNormal(double x, double mu, double sigma)
    {
	return 0.5 + 0.5 * erf((x - mu)/(sigma * Math.sqrt(2.0)));
    }

    // erf function taken from:
    // http://introcs.cs.princeton.edu/java/21function/ErrorFunction.java.html
    // fractional error in math formula less than 1.2 * 10 ^ -7.
    // although subject to catastrophic cancellation when z in very close to 0
    // from Chebyshev fitting formula for erf(z) from Numerical Recipes, 6.2
    public static double erf(double z) {
        double t = 1.0 / (1.0 + 0.5 * Math.abs(z));
        // use Horner's method
        double ans = 1 - t * Math.exp( -z*z   -   1.26551223 +
				       t * ( 1.00002368 +
					     t * ( 0.37409196 + 
						   t * ( 0.09678418 + 
							 t * (-0.18628806 + 
							      t * ( 0.27886807 + 
								    t * (-1.13520398 + 
									 t * ( 1.48851587 + 
									       t * (-0.82215223 + 
										    t * ( 0.17087277))))))))));
        if (z >= 0) return  ans;
        else        return -ans;
    }


}